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Abstract. We address the problem of verifying timed properties of 
Markovian models of large populations of interacting agents, modelled 
as finite state automata. In particular, we focus on time-bounded prop¬ 
erties of (random) individual agents specified by Deterministic Timed 
Automata (DTA) endowed with a single clock. Exploiting ideas from 
fluid approximation, we estimate the satisfaction probability of the DTA 
properties by reducing it to the computation of the transient probability 
of a subclass of Time-Inhomogeneous Markov Renewal Processes with 
exponentially and deterministically-timed transitions, and a small state 
space. For this subclass of models, we show how to derive a set of De¬ 
lay Differential Equations (DDE), whose numerical solution provides a 
fast and accurate estimate of the satisfaction probability. In the paper, 
we also prove the asymptotic convergence of the approach, and exem¬ 
plify the method on a simple epidemic spreading model. Finally, we also 
show how to construct a system of DDEs to efficiently approximate the 
average number of agents that satisfy the DTA specification. 

Keywords: Stochastic Model Checking, Fluid Model Checking, Deter¬ 
ministic Timed Automata, Time-Inhomogeneous Markov Renewal Pro¬ 
cesses, Fluid Approximation, Delay Differential Equations. 


1 Introduction 

One of the major technological challenges in computer science and engineering 
is the design and analysis of large-scale distributed systems, where many au¬ 
tonomous components interact in an open environment. Examples include the 
public and shared transportation in smart cities, the power distribution in smart 
grids, and the robust communication protocols of online multimedia services. 
In this context, the mathematical and computational modelling plays a crucial 
role in the management of such Collective Adaptive Systems (CAS), due to the 
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need of understanding and control of their emergent behaviours in open work¬ 
ing conditions. The intrinsic uncertainty of CAS can be properly captured by 
stochastic models , but the large number of interacting entities always results in a 
severe state space explosion , introducing exceptional computational challenges. 
In particular, the scalability of the models and of their analysis techniques is 
a major issue in the development of stochastic model checking procedures for 
the verification of formal properties. In this context, up to now, the numerical 
approaches 24 are deeply hampered by the state space explosion of the large 


stochastic models, and the statistical methods based on simulation require a 
large computational effort. 

A recent line of work tries to address the issue of scalability by exploit¬ 
ing stochastic approximation techniques 10 11 , like the Fluid Approximation 
J8jj9j[l8|. In this method, a stochastic discrete model is replaced by a simpler 
continuous one, whose dynamics is described by a set of differential equations. 
In [8], the authors exploit this limit construction to verify properties that asses 
the behaviour of a single individual in a collective system, and define a procedure 
called the Fluid Model Checking (FMC) 17 25 . This technique is based on the 
Fast Simulation Theorem 161, which ensures that in a large population, a single 
entity is influenced only by the mean behaviour of the rest of the agents. 

In this work, we extend | 8] to more complex time-bounded properties specified 
by Deterministic Timed Automata endowed with a single clock l,|3j 17 . As 
in mmm- we combine the agent and the DTA specification with a product 
construction, obtaining a Time-Inhomogeneous Markov Renewal Process |15| . 
We then exploit results [6,22, defining the Fluid Approximation of this type of 
models as the solution of a system of Delay Differential Equations (DDE ) |16| . 
Other works dealing with the verification of DTA properties are [4|[l2j[n|[T9p 

Main Result. We introduce a new fast and efficient Fluid Model Checking 
procedure to accurately approximate the probability that a single agent satisfies 
a single-clock DTA specification up to time T. Similarly to [8], the technique is 
based of the Fast Simulation Theorem, and couples the Fluid Approximation of 
the collective system with a set of Delay Differential Equations for the transient 
probability of the Time-Inhomogeneous Markov Renewal Process obtained by 
the product construction between the single agent and the DTA specification. 

In the paper, we discuss the theoretical aspects of our approach, proving 
the convergence of the estimated probability to the true one in the limit of an 
infinite population. We also show the procedure at work on a running example 
of a simple epidemic process, emphasising the quality of the approximation and 
the gain in terms of computational time. Finally, by exploiting the construction 
of 10 22 , we also show how to define a set of DDEs approximating the mean 


number of agents satisfying a single-clock DTA specification up to time T. 

Paper structure. In Sec. [2j we introduce the modelling language, the Fluid Ap¬ 
proximation, the Fast Simulation Theorem, and the DTA specification for the 
timed properties. In Sec. [3j we present our FMC procedure, defining the DDEs 


for the probability that the single agent satisfies the timed property. In Sec. 3.1 


we adapt our verification technique to compute the mean number of agents that 
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meet the DTA requirement. In Sec. [4j we discuss the quality of the approxima¬ 
tion on the epidemic example. Finally, in Sec. [5] we draw the final conclusions. 
The proofs of the theoretical results are reported in the Appendix. 


2 Background and Modelling Language 


Agent Classes and Population Models. A collective system is comprised 
of a large number of interacting agents. To describe its dynamics, we define a 


population model 


in which the agents are divided into classes, called agent 


classes, according to their behaviour. 


Definition 1 (Agent Class). An agent class srf is a pair ( S , E ) in which S = 
{1,..., m} is the (finite) state space and E = {ei ,..., e v } C S x Jz? x S is the 
(finite) set of local transitions of the form = Sj s[, where Si,s[ £ S are the 
initial and arrival states, and ai £ Jz? is the unique label of a, i.e. a* A otj for 

i ^£1 


Intuitively, an agent in class = ( S , E) is a finite state automaton that 
can change state by performing the actions in E. Then, assuming that agents 
in the same state are indistinguishable, to define the population model, we rely 
on the counting abstraction, counting how many agents are in each state at time 
t. Hence, for each agent class, we define the collective or counting variables 
x[ N \t),..., Xm\t) given by X i j N \t) = l {F w (t)=j} , where Y k (N) (t) £ 

{1,..., m} is the random variable denoting the state of agent k at time t, and 
N = y irf yr XJ N ^ is the population size, that we assume constant in time (cf. 
also jl]). Then, given n = y^|5|, the state of the population model is given 
by the vector X' A b (t) £ (It>o) n that enlists the counting variables of the agent 
classes. 


Definition 2 (Population Model). A population model X < - N ' 1 is a tuple = 
(A, T (n \x{, N) ), where A = {sf \,..., srf v } is the set of agent classes, as in Defi¬ 
nition 0 = X^)(0) is the initial state; and = {n,... ,iy} is the set 

of global transitions of the form Ti = (Sj, f\ N \ v,[ JV ' ) ), where: 


— Si = {|si —A s[,... ,s p —A Sp|} is the (finite) multi-set of local transitions 
synchronized by Tt; 

— f*> N ^ : (R> 0 )" —> lR>o is the (Lipschitz continuous) global rate function; 

— Vi = Sa es-l{l a il}l(-®-s3 — lg'.) is the update vector, where |{|oj|}| is the mul¬ 
tiplicity of oti in Si, and l Sj . is the vector equal to 1 on Sj and 0 elsewhere. 


The restriction on the uniqueness of the labels can be dropped (as in 10 ) at the 


price of heavier notation and combinatorics in the definitions of the rest of the paper. 
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Fig. 1. The agent class si (left) and property D (right) of the running example. 


When a global transition Ti = (§,;, f N \ Vj) is taken, the transitions in Si fire at 
the local level, meaning that, for each s s' in Si, an agent moves from s to 
s'. Hence, the update vector v* encodes the net change in the state X (A b(f) of 
X l ' N> due to transition r,. Moreover, for the model to be meaningful, whenever 
at time t it is not possible to execute t,, because there are not enough agents 
available, i.e. (X^(t) — v») < 0 for some j £ {1,... ,n} with n = |X^(t) |, 

we require the rate function to be zero, i.e. fl N \~K^ N ^(t)) = 0. 

Example. The running example that we consider is a simple SIS model, describ¬ 
ing the spreading of a disease inside a population. All agents belong to the same 
agent class si, depicted in Fig. |TJ and can be either susceptible (S) or infected 
(I). When they are susceptible , they can be infected (inf), and when they are in¬ 
fected, they can either pass the infection (pass) or recover (rec). Hence, the state 
X^ (t) of the population model is XW (t) = (Xg N ' > (t), Xj N) (t)), and we define 2 

global transitions: T r = ({I AES. f N ^) and r, = ({S -^4 I, I pass > f( N 1). 
The former, r r , mimics the recovery of one entity inside the population, while 

Ti synchronises two local actions, namely S —A I and I ——> I, and mod¬ 
els the transmission of the virus from an infected agent to a susceptible one. 
Finally, the rate functions depend on the number of agents involved in the tran¬ 
sitions and follow the classical rule of mass action fr N \t) = k r Xj N \t) and 
fi N \t) = jfkiXg N \t)Xj N \t), where k r ,ki £ K> 0 . 


Fluid Approximation. The Fluid Approximation [s}[9 18 of a population 
model AW = (A ,T {N) ,x. { Q N) ) is an estimate of the mean behaviour of its agents. 
To compute this approximation, we first normalise X ( iV * by dividing the state 
vector XS N \t) and the initial state xj^ by the population size N, i.e. we define 
XW^) = ~xS N \t)/N and x < ' 0 N ' 1 = x' 0 N ' 1 /N. Then, for all transition £ T^ N \ 
we let fl^ (X) be the rate function, where we substitute the counting variables 
of X. iN) (t) with the new normalised counting variables of X(i). Moreover, we 
assume that for each f N ' l (X.), there exist a Lipschitz function fi(X.) such that 
> /j(X) uniformly. Finally, in terms of _/)(X), we define the 


r\x)/N - 


-+oo 


drift F(X) given by F(X) = ]T r . v it fj(X), whose components represent the 
instantaneous net flux of agents in each state of the model. Then, given a time 
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horizon T < +oo, the Fluid Approximation t ) of X^ N ^ is the uniqu^] solution 
of the system of Ordinary Differential Equations (ODEs) given by 




for 0 < t < T, 


with #(0) = xo- The accuracy of the approximation improves the larger is the 
ensemble of agents that we consider, i.e. the larger is TV, and is exact in the limit 
of an infinite population. Indeed, the following theorem holds true 
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Theorem 1 (Fluid Approximation). For any T < +oo and e > 0, 


Prob 


sup ||xW(t)-#(t)||> 
0 <t<T 


N->-\-oo 


> o. 


Fast Simulation. In this paper, we are interested in the behaviour of a (ran¬ 
dom) single agent inside a population. As we have just seen, the dynamics of 
a large population can be accurately described by a deterministic limit, the 
Fluid Approximation. But when we focus on one single agent in a collective sys¬ 
tem, we need to keep in mind that its behaviour in time will always remain a 
stochastic process, even in large populations. Nevertheless, the Fast Simulation 
Theorem jh, 16 20| guarantees that in the limit of an infinite population size, 
the stochastic process of the single agent senses only the mean behaviour of the 
rest of the agents (i.e. there is no need to keep track of all the states of all the 
other entities in the population). This means that, when the population size is 
large enough, to analyse the dynamics the single agent, we can define the Fluid 
Approximation of the population model, and then use its state (i.e. the mean 
state of the rest of the agents) to compute the rates of a Time-Inhomogeneous 
CTMC (ICTMC) pi that describes the behaviour of the single agent. 

Formally, let F™) (t) be the stochastic process that describes the state of the 
single agent in the population model X^ = (A, , Xq^ ) with state vector 

XW(t). By definition, Y^ N \t) is not independent of Now consider 

the normalised model X ^ described by XS N \t), and let <P(t) be the Fluid 
Approximation of X^ N \ Define the generator matrix Q {N) (x) = {q { ip(x)) of 
Y( N \t) as a function of the normalised counting variables, i.e. V q[^\x), 

Prob {F (JV) (t + dt) = j | Y^ N \t) = i, X^ft) = x} = q\f\x)dt. 

Notice that Q (N \x) can be computed from the rates in X^ N \ Indeed, for i j, 


dij 


(*) = 

t€T l 


Xj 


N 


where |{|z —>• j € § T [}| is the multiplicity of i —> j in the transition set § r of r, 
i.e. the number of agents that take such transition according to r. Furthermore, 

4 Existence and uniqueness of d>(t) are guaranteed by the Lipschitzianity of the /i(X). 
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as customary, let q\f\x) = ')■ Then, since f} N) (X)/N N ^ f+ T > > 

/i(X), we have that Q l ' N \x) —> Q(x), where Q{x) is computed in terms of the 
Lipschitz limits _/)(X). Now, define the stochastic processes: 

1. ZW(t), that describes the state of the process Y^ N \t) for the single agent 
in class sY, when Y^ N \t) is marginalised from 

2. Z(t), that is the ICTMC, defined on the same state space of Z^ N \t), with 
time-dependent generator matrix Q(<P(t)), i.e. the generator matrix Q(t), 
where the Lipschitz limits fi(t) are computed over the components of <P(f). 


Then, the following theorem can be proved 
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Theorem 2 (Fast Simulation). For any time horizon T < +oo and e > 0, 


Probl sup \\zW(t)-Z(t)\\> 

1.0 <t<T 


N—t+oo 


> 0 . 


Example. For the running example, if we consider a population of 1000 agents, i.e 
N = 1000, and an initial state = (900,100), then the Fluid Approximation 
<P(t) of the population model is the unique solution of the following ODEs: 

i^(t) = -k i ^i(t)^s(t) + k r $ I (ty, (<P S ( 0) = 0.9; 

1 ! %{t) = +k i $ I {t)$s{t)-k r $ I (ty, \«p/(o) = o.i. 

The generator Q(4>(t)) of the ICTMC Z(t) for the single agent, instead, is: 
qs,s{t) = ~q S ,i(ty, qsj(t) =k i $ I (t); q I ,s{t)=k r ; q ItI (t) =-q I>s (t). ( 2) 


2.1 Timed Properties 

We are interested in properties specifying how a single agent behaves in time. 
In order to monitor such requirements, we assign to it a unique personal clock , 
which starts at time 0 and can be reset whenever the agent undergoes specific 
transitions. In this way, the properties that we consider can be specified by a 
single-clock Deterministic Timed Automata (DTA) 0[f3], which keeps track of 
the behaviour of the single agent with respect to its personal clock. Moreover, 
since we want to exploit the Fast Simulation Theorem, we restrict ourselves to 
time bounded properties and, hence, we assign to the DTA a finite time horizon 
T < Too, within which the requirement must be true. 

Definition 3 (Timed Properties). A timed property for a single agent in 
agent class sY is specified as a single-clock DTA of the form D = O(T) = 
(T, «Sf, c,CC, Q,qo, F,—t), where T < Too is the finite time horizon; JZ is the 
label set of sY; c is the personal clock; CC is the set of clock constraints, which 
are conjunctions of atoms of the form c < X, c < A, c > X or c > X for A € Q; Q 
is the (finite) set of states; qo £ Q is the initial state; F C Q is the set of final 
(or accepting) states; and — > C Q x x CC x {0,{c}} x Q is the edge relation. 
Moreover, D has to satisfy: 
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— (determinism) for each initial state q £ Q, label a £ Jz?, clock constraint 
Cpx, £ CC, and clock valuation rj{c) £ R>o, there exists exactly one edge 
q ’ - r > q' such that r]{c) \=cc Ct<^}' 

— (absorption) the final states are all absorbing. 

A timed property D is assessed over the time-bounded paths (of total duration 
T) of the agent class si sampled from the stochastic processes Z^ N \t) and Z{t) 
defined for the Fast Simulation in Sec. [2j The labels of the transitions of si 
act as inputs for the DTA O, and the latter is defined in such a way that it 
accepts a time-bounded path a if and only if the behaviour of the single agent 
encoded in a satisfies the property represented by D. Formally, a time-bounded 
path a = s 0 Si ... a ’ 1 ’* - > s n+ \ of si sampled from Z^ N \t) (resp. 

Z{t)), with £"= 0 tj < T, is accepted by D if and only if there exists a path 
go -^4 gl 1 ) —4 g( 2 ) —4 ... -^4 q( n + 1 ) of O such that q( n+1 l £ F. In the path 
of D, gb +1 ) £ Q denotes the (unique) state that can be reached form gb) £ Q 
taking the action gb) gb+L whose clock constraint is satisfied by 

the clock valuation 77 (c) updated according to time ti. In the following, we will 
denote by ,t the set of time-bounded paths of si accepted by D. 

Example. We consider the following property for the running example: within 
time T, the agent gets infected at least once during the A = 5 time units that 
follow a recovery. To verify such requirement, we use the DTA D = O(T) rep¬ 
resented in Fig. [T| If we record the actions of the single agent on D, i.e. we 
synchronise si and D, when the agent recovers (rec), D passes from state gO to 
gl, resetting the personal clock c. After that, if the agent gets infected {inf) 
within 5 time units, the property is satisfied, and D passes from state gl to g2, 
which is accepting. If instead the agent is infected {inf) after 5 units of time, 
D moves back to state gO, and we start monitoring the behaviour of the agent 
again. In red we highlight the transition that resets the personal clock c in D. 

3 Fluid Model Checking of Timed Properties 

Consider a single agent of class si = {S,E) in a population model = 

{si,T( N>> jXq^), and a timed property D = O(T) = (T, i£, Ig, CC, Q, go, F, —►). 
Let Esj ,®,t be the set of time-bounded paths of si accepted by D. Moreover, let 
ZW(i) and Z(t) be the two stochastic processes defined for the Fast Simulation 
in Sec. [2j The following result holds true. 

Proposition 1. The set AUqn ,t is measurable for the probability measures 
Prob Z (N) and Probz defined over the paths of Z^ N \t) andZ{t), respectively. □ 

Let pW)(T) = Prob Z ( jv){X'^ ) b,t} and P{T) = Probz{£srf,o,T}- In this paper, 
we are interested in the satisfaction probability p( N l{T), i.e. the probability 

5 The notation 77 (c) |=cc Cx stands for the fact that the value of the valuation 77 (c) of 
c satisfies the clock constraint c M . 
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that the single agent satisfies property D within time T in X^ N l . Then, the main 
result that we exploit in our Fluid Model Checking procedure is that, when the 
population is large enough (i.e N is large enough), pW(T) can be accurately 
approximated by P(T), which is computed over the ICTMC Z(t), whose rates 
are defined in terms of the Fluid Approximation &(t) of X^ N \ The correctness 
of the approximation relies on the Fast Simulation Theorem and is guaranteed 
by the following result. 

Theorem 3. For any T < +oo, lim^v-^oo P^(T) = P(T). □ 

Moreover, to compute P{T ), we consider a suitable product construction = 
srf ® D, whose state is described by a Time-Inhomogeneous Markov Renewal 
Process (IMRP) [l5] that we denote by In the rest of this section, we 

define and Z^ n {t), and we show how to compute the satisfaction probability 
P(T) in terms of the transient probability P(T ) of Z s / a {t'). 

The Product s$q. We now introduce the product between srf and D, whose 
state is described by a Time-Inhomogeneous Markov Renewal Process (IMRP) 
Z,j K (t) that has rates computed over the Fluid Approximation &(t) of X^ N \ 

A Markov Renewal Process (MRP) [15] is a jump-process, where the sojourn 
times in the states can have a general probability distribution. In particular, in 
the MRP Z^ft), we will allow both exponentially and deterministically-timed 
transitions, and in the following, we will refer to them as the Markovian and 
deterministic transitions , respectively. Since the transition rates of Z^ n {t) will 
be time-dependent, Z^ B (t) will be a Time-Inhomogeneous MRP. 

To define the product So, {-M, £}, so,b, ^b)) let <5i < ... < 6k be 

the (ordered) constants that appear in the clock constraints of D, and extend 
the sequence with <5o = 0 and = T. The state space So of s$o is given by 
{ 1 , ..., fc + 1} x S x Q. The first element of 5 b identifies a time region of the 
clock c, and we refer to 5o* = {(*, s, q) \ s G 5, q G Q} as the i-th Time Region 
of 5 b- The rest of s$o will be defined in such a way that the agent is in Sn i if 
and only if c satisfies 8i_\ < 77 (c) < 5i, where r] is the valuation of c. 

The set A4 of Markovian transitions of is the smallest relation such that 


V i 1,..., k T 1 , 


s' £ E A q > q’ A [£*_!,&] \= 

(i, s, q) A ( i,s',q') G M 


(3) 


V i G 1,..., k H - lj 


s' G E A q --- > q' G —> A [di_i,di] |= Cxi 

( i,s,q ) A (1 ,s',q') G M 


(4) 


Intuitively, rule (|3j) synchronises the local transitions s A- s' G E of the 
agent class = ( S,E ) with the transition q a ’ c c<1 ’ 0 > q' G— >• that has the same 
label in D, obtaining a local transition ( i,s,q ) A ( i,s',q') G M in for each 
time region i whose time interval [Si-i,6i] C [0, T] satisfies the clock constraint 
ex,, meaning that Vf G [8i- i,<5i], t \= Cx- Rule Q, instead, defines the reset 
transitions (i, s, q) —> (1, s' , q') G A4 that reset the personal clock c either within 







the I s * Time Region (when i = 1), or by bringing the agent back to the 1 st Time 
Region. In the following, we denote by 1Z C M. the set of the reset transitions. 

To describe the deterministic transitions of g/g, instead, we define a set £ of 
clock events. Each clock event has the form e = (A e , A,p e ), where A e C Sb ; is 
the active set , A is the duration , and p e : A e x So —► [0,1] is the probability 
distribution. If the agent enters A e , that is the sets of states in which e can be 
active, a countdown starts from A. When this elapses, e, is deactivated and the 
agent is immediately moved to a new state sampled from p e ((i , s, q), •) : Sjj —> 
[0,1], where (z, s, q) £ A e is the state in which the agent is when the countdown 
hits zero. Moreover, e* is deactivated also when the agent takes a reset transition. 
In ^4), we define: 

— one clock event e, £ £ for each time region Sm> * = 2,..., k; 

— £ + 1 clock events e°, e \,..., e\ £ £ for the 1 st Time Region, where £ is the 
number of reset events (1, s, q) (1, s', q') £ 1Z defined by (Jr]) with z = 1. 

For i > 1, Ai = Sm , Ai = Si — 4-1, and the probability distribution is 


Pi{{i,s,q ), ( i',s\q')) 


1 if i! = i + 1, s' = s,q' = q 1 
0 otherwise. 


(5) 


As it is defined, each clock event e* with i > 1 connects each state (z, s, q) £ A, 
with (z + l,s, g) £ Sbi+i, hence, when the duration A,; of e* elapses, the clock 
event moves the agent from its state to the equivalent one in the next time 
region. When z = 1, instead, the duration and the probability distribution of 
each clock event e{, j = are defined in the same way as before (i.e. 

Aj = <5i — <5o = <5i, and p{ is given by ([H])) , but the activation sets are now 
subsets of Sb!• Indeed, since each reset transition (l,s,g) —4 (1 ,s',q') £ 1Z 
initiates the clock, for each of them, we need to define a clock event ef, whose 
activation set Aj is the set of states in 5 bj that can be reached by the agent 
after it has taken the reset transition (l,s,g) —A (1, s', <j ,/ ). Furthermore, we 
have to define an extra clock event e°, with A\ = Sbi, A^ = S -\, and p\ given 
by 0, that is the only clock event initiated at time t = 0 (and not by the 
agent entering *4°). Indeed, we require for the initial state so,b of s$o to be 
one of the states of the form (l,s,g 0 )j where s £ S and qo is the initial state 
of D (hence, So,b belongs to A\). Finally, since the probability distributions 
p {, Vj, are all defined as in also the clock events of the 1 st Time Region 
move the agent from a state to the equivalent one in the next time region (the 
2 nd ), when the countdown from A\ = <5i elapses. In the following, we denote 
by (i,s,q) ~—* e (z + 1 ,s,q) the deterministic transition from ( i,s,q ) £ Sm to 
(z + 1 ,s,q) £ S]m + i encoded by e £ £, and by v e ,s,q = l(i+i,«, 9 ) - its 

update vector. The last component of gfa that we define is the set of final states 
Fji> which is given by F a = {(z, s,q) £ 5b | q £ F}. 


Example. Fig. [2] represents the product srf \j between the agent class and the 
property D of the running example (Fig. [l]). The state (1,7, ql) that cannot be 
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1 st Time Region 
0 < c< 5 


2 nd Time Region 
c > 5 



Fig. 2. The agent class s/o associated with the DTA B of the running example. 


reached by the single agent is omitted. The black transitions are the Markovian 
transitions without reset; the red transitions are the Markovian transitions that 
reset the clock; finally, we define 2 clock events, e® and ej , with duration A = 5 
for the I s * Time Region, and the dashed green (resp. blue) transitions are the 
deterministic transitions encoded by e° (resp. e\). In blue, we also highlight the 
states that belong to the activation set A e i (while A e o is the whole 1 st Time 
Region). Intuitively, the agent can be found in one of the states belonging to the 
1 st Time Region whenever its personal clock c is between 0 and 5, i.e. less that 5 
time units have passed since t = 0 or since a recovery rec. In a similar way, the 
agent is in the 2 nd Time Region when the valuation of c is above 5. Moreover, 
when the the duration of the clock events elapses (i.e. the countdown from 5 hits 
0), the agent is moved from the 1 st Time Region to the 2 nd Time Region by the 
deterministic green and blue transitions, that indeed have duration A = 5 and 
are initiated at t = 0 or by the reset actions rec, respectively. At the end, the 
agent is in one of the final states ((1,5, q2), (1, Z, ^2), ( 2 , 5 , 92 ) or (2, 1, q2)) at 
time T, if it meets property D within time T, i.e. within T, the agent has been 
infected during the 5 time units that follow a recovery. Hence, to verify D, we 
will compute the probability of being in one of the final states of s/d at time T. 


The IMRP Z#? D (t) and the Satisfaction Probability P(T). Now we show 
how to formally define the IMRP Z^ 0 (t) that describes the state of the product 
s/d in the mean field regime. In particular, we derive the Delay Differential 
Equations (DDE) |l5] for the transient probability Pit) of Z^(t), in terms of 
which we compute the satisfaction probability P{T). 


Let <P(f) be the Fluid Approximation of the population model . To define 

the transient probability P(f) of Z^ B {t), we exploit the fact that, in the case 
of an IMRP, we have: ^£(i) = M{&{t))P{t) + D($(t),P(t)) (cf. 15 ). In this 


equation, M(4>(t)) is the generator matrix for the Markovian transitions, and 
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D(<P(t), P(t)) accounts for the deterministic events. The elements of 
are computed following the same procedure that was described in Sec. [2j where 
the multiplicity of each transition (i, s, q) ( i , s', q') £ Ai in is always equal 
to 1 (one single agent) and the Lipschitz limit / Q (#(£)) of a is that of the rate 
of the transition s ^ s' in X^ N ^ from which a was derived (by rules ^ or §)• 

To define the components of D(&(t), P(t)), instead, consider any clock event 
e = (Ai, Ai,pi) £ £, except e®, whose contribute will be computed later oi 0 
Choose one of the deterministic transitions (i, s , q) ~—* ei (i + 1, s, q) encoded by 
e*. The agent takes this transition at time t when: (1) it entered A t C So i at 
time t — Ai (initiating its personal clock), and (2) it is in state ( i,s,q ) £ Ai 
at time t (when the duration of ej elapses). Hence, to compute the term that 
corresponds to this transition in £)(#(£), P(t)), we need to: (1) record the flux of 
probability that entered Ai at time t— Ai, and (2) multiply it by the probability 
that the agent reaches ( i,s,q ) £ Ai at time t, conditional on the state at which 
it entered A, at t — Ai. 

To compute the probability of step (2), we need to keep track of the dynamics 
of the agent while the clock event e,; is active. For this purpose, let Ai be the 
activation set Ai of e, extended to contain an extra state s out = (i, s out ,q ou t), 
and let Ai be the set Ai of Markovian transitions in £/o modified in order 
to make the reset transitions that start in Ai finish into s ou t (i.e. for every 
(i , s, q) A- ( i', s ', q') £ 1Z C Ai, we define (i, s, q) A- s ou t £ Ai), and to have s ou t 
absorbing Let G i(<P(t)) £ Matr(|*4j| x |^4j|) be the time-dependent matrix s.t. 


(G i(^(0))(i, s ,g),( M ', ? ') - 


i,s,q)- 


E 

->■(&,s',g')E.M 


Mt) 


urn) 


( 6 ) 


where again the Lipschitz limit f a (t) of each a £ Ai is that of the transition s —> 
s' in from which its copy a £ Ai was derived (by ^ and 0 ) . Moreover, 

let the diagonal elements of G i(<P(t)) to be defined so that the rows sum up 
to zero. Then, we introduce the probability matrix Yi(t ), which is computed in 
terms of the generator G* (#(£)) according to the following ODEs (see also 

f(t) = Yi(t)Gi(m) - Gm - Ai))Yi(t), Ai<t< T, 

\^(t) = Yi(t)Gi($(t)), 0 <t<Ai, 

with 'Kj(O) = I. By definition, {Yi{i))u^ s i^ is the Fluid Approximation of 

the probability of step (2), i.e. the probability that the agent, which has entered 
Ai in state (*, s' , q') at time t — Ai, moves (Markovianly) within Ai for Ai units 
of time, and reaches (i,s,q) £ Ai at time t (exactly when e, elapses). 

In terms of the probability matrix Yi(t), we can now define the component 
of P(t)) that corresponds to the deterministic transition (i,s,q) —-> ei 

6 If e is one of events of the 1 st Time Region, i.e. e = e(, for some j — 1, ...,£, in this 
section, we drop the index j to ease the notation, i.e. we write e\ = e\ = (Ai, Ai,pi). 

1 The absorbing state s ou t is needed for the probability Yi ( t ) of step (2) to be well 
defined. Indeed, the agent can deactivate e, by taking a reset transition. 
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(i + 1, s, q) of the clock event e, £ £. This component is the element in posi¬ 
tion ((z, s, q), ('i + 1, s, q)) in P{t)), we call it D ei>Sjq (&(t), P(t)), and is 

given by 


D euS ^{t),P{t))= y 

(i,S,q)eAi 


!{»>! }D ei _ u s,q{&{t - Ai),P{t 


Ai)) + l{i=l} x 


E 

(i',s' ,q') — y(l,s,q)£lZ 






(^i(i))(i,s, q),(i,s.q)t 


( 8 ) 

where (P(t — A\ ))(i^ s ',q') is the component in position {i' ,s' ,q') £ So., in the 
vector of the transient probability P(t — Ai) of Z^ B at time t — A\. In ([8]), 
for each state (z,s, q) in the activation set Ai, the quantity inside the squared 
brackets is the probability flux that entered (z, s, q) at time t — Ai. In particular, 
when i > 1, D ei _ lt s,q{&{t — Ai),P(t — Ai)) accounts for the termination of 
clock event e^_i (i.e. the deterministic transition (z — 1 ,s,q) —-» e< ( i,s,q ) fired 
at time t — Ai). When we consider the 1 st Time Region, i.e. i = 1, instead, 
each term in the sum over the reset transitions is the flux of probability entering 
(1 ,s,q) at time t — A\ due to a clock reset. Finally, (T"j(t))(i,s,g),(i,s, 9 ) again 
the probability of reaching (z, s, q) £ Ai from (z, s, q) £ Ai in Ai units of time. 

All the other off-diagonal elements of D(4>(t), P(t)) can be computed in a 
similar way, while the diagonal ones are defined so that the rows sum up to zero. 
Moreover, since at the end D(&(t),P(t)) depends on the state of the system 
at times t — A 1 ,... ,t — A k (through the probabilities Yi(t), i = 1 ,..., k), we 
write D(<P(t)) = D(<1>, P,A\,... , Ak, t). Then, we define the transient probabil¬ 
ity P{t) of the IMRP Z^it) as the solution of the following system of DDEs: 


P(t) = f M(s)P(t)ds +[ D(&,P,A 1 ,...,A k ,s)ds + l t > Al Y Ve^e^a.r 

(s,q)eSxQ 

(9) 


In@, the third term is a deterministic jump in the value of Pit) at time t = A\, 
and represents the contribute of the clock event ej. In such term, the vectors 
i'e°,s,q are the update vectors of the deterministic transitions encoded by ej 
(hence, the sum is computed over all such transitions), and the probability y e o 
is the value at time t = A\ of the component in position (so,b, (1, s, q)) (where 
sq,!} is the initial state of &/o) in the matrix Y e o(t) defined by: 


dY e 

dt 


Ht) = Y e a(t)G 1 (#{t)), 


0<t<A u 


with Gi(^(t)) defined as in and y e o(0) = I. Hence, 2/ e ;=C^e;(^i)) S o,„,(i,s,<r) 
is the probability that, starting form so,b, the agents reaches (l,s, <?) £ Sox at 
time t = Ai (exactly when the deterministic event (1 ,s,q) ---> e o (2, s,q) fires). 

Given the product s^o, the IMRP Z^ B (t), and its transient probability P(t), 
the following result holds true. 
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Proposition 2 . There is a 1:1 correspondence between E^ t o,T and the set 
AccPath(srfo,T) of accepted paths of duration T of srfo- Hence, 

P(T) = Prob z {E^,o, T } = Prob Zs j n {AccPath{srf 0 ,T)} = Pf b (T), 

where Probz^ B is the probability measure defined by Z^ D , and Pf„{T) is the sum 
of the components of P(T) corresponding to the final states Pd of s&q- □ 

In other words, according to Proposition^ when the population of X W is large 
enough, Pp a (T) is an accurate approximation of the probability that a (random) 
single agent in X satisfies property D within time T. 

Example. For the product s^b i n Fig- [2j the non-zero off-diagonal entries of 
the generator matrix G e i(&(t)) of the clock event e\ are: G^s, q i)(i, q 2 )(t) = 
ki$i(t)\ G(s,q 2 )(i,q 2 )(t) = ki$i(t )\and G(j, g2 )(s,g 2 )0) = K. In terms of G e i 
we can define Y e i(t), as in (Pfj) , that is then used in the DDEs d9l) for the prob¬ 
ability P(t). In this latter sex of 9 DDEs (one for each state of £%), we have: 


P(i,s, q i){t) = / k r P(i,s, q i)(s)ds - kid> I (s)P( lt s,qi){s)ds + 

Jo Jo 

— / k r Y(l'S,ql),(l,S,ql)( S — 5, s)P(l } S : ql)(s)dS. 

Jo 

Remark. The presence of only one clock in D enables us to define srfo in such a way 
that Zjz/ B (t) is an IMRP. This cannot be done when we consider multiple clocks 
in D. Indeed, in the latter case, the definition of the stochastic process which 
describes the state of the product .s/n is much more complicated, since, when a 
reset event occurs, we still need to keep track of the valuations of all the other 
clocks in the model (hence, the dynamics between the time regions of srf \j is not 
as simple as in the case of one single clock). In the future, we plan to investigate 
possible extensions of our model checking procedure to timed properties with 
multiple clocks, also taking into account the results of 19 and [4|. 


3.1 The Mean Behaviour of the Population Model 

It is possible to modify our FMC procedure in order to compute the mean num¬ 
ber of agents that satisfy D. This can be done by assigning a personal clock 
to each agent, and monitoring all of them using as agent class the product s^o 
defined in Sec. [3j In terms of ^d, we can build the population model Xg, with 
s^b as the only agent class, and the sum Pp B (T) of the components correspond¬ 
ing to the final states of srfa in the Fluid Approximation &(t) of Xo computed 
at t = T is indeed the mean number of agents satisfying property D within 
time T. The construction of Xb is not difficult: it follows the procedure of 110 , 
where a little extra care has to be taken just in the definition of the global 
transitions of Xb- Indeed, if we build for instance the population model Xo 
of the running example, we need to consider that the infected individual that 
passes the virus to an agent in state (l,S,qO) can be now in one of five states: 
(1, 1, qO), (1, /, q2), (2, J, qO), (2,1, ql) or (2,/,q2). For this reason, we have to 


13 



Fluid Model Checking 


N 

MeanRelErr 

MaxRelErr 

RelErr(T) 

TimeDES 

TimeFMC 

Speedup 

250 

0.0927 

6.4512 

0.1043 

11.0273 

0.4731 

23.3086 

500 

0.0204 

1.7191 

0.0048 

44.0631 

0.3980 

110.7113 

1000 

0.0118 

0.7846 

0.0003 

170.9154 

0.3998 

427.5022 


Fluid Approximation of the mean behaviour 


N 

MeanRelErr 

MaxRelErr 

RelErr(T) 

TimeDES 

TimeFluid 

Speedup 

250 

0.1127 

0.2316 

0.0921 

105.5647 

0.4432 

339.7217 

500 

0.0289 

0.3177 

0.0289 

415.0635 

0.4237 

979.6165 

1000 

0.0117 

0.2216 

0.0117 

1547.0340 

0.4213 

3672.0484 


Table 1. Mean Relative Error (MeanRelErr), Maximum Relative Error (MaxRelErr), 
and Relative Error at final time (RelErr(T)) of the FMC (top) and the Fluid Approx¬ 
imation of the mean behaviour (bottom) for different values of N. The table enlists 
also the execution times (in seconds) of the DES (TimeDES) and the approximations 
(TimeFMC and TimeFluid), and the speedups (TimeDES divided by the other times). 


define five Markovian global transition in Xjj, each of which moves an agent 
from (1,5, qO) to (l,/,qO) at a rate that is influenced by the number of indi¬ 
viduals that are in the infected states of recorded in the counting variables 
X(i t i iq o)(t),X( ltIiq 2 )(t),X( 2 ,i,qO)(t),X( 2 ,i,qi)(t) or X( 2 ,i,q 2 )(t). The same reason¬ 
ing has to be followed for the definition of the infections of the agents in states 
(1, 5, ql), (1, 5, q2), (2, 5, qO), (2, 5, q\) and (2, 5, q2). At the end, as for the single 
agent, due to the deterministic events, the Fluid Approximation 4>(t) of Xo is 
the solution of a system of DDEs similar to Q. The definition of such approx¬ 
imating equations for a population model with exponential and deterministic 
transitions is not new 22 , but, even if the results are promising (see Sec. El), to 


our knowledge, nobody has yet proven the convergence of the estimation in the 
limit N —> +oo. We save the investigation of this result for future work. 


4 Experimental Results 

To validate the procedures of Sec. |3j we performed a set of experiments on 
the running example, where we fixed: ki = 1.2, k r = 1, A = 5, and an initial 
state of the population model with a susceptible-infected ratio of 9:1. As in 
Fig-! we let the single agent start in the susceptible state, and we considered 
three different values of the population size: N = 250, 500,1000. For each N, we 
compared our procedures with a statistical estimate from 10000 runs, obtained 
by a dedicated Java implementation of a Discrete Event Simulator (DES). The 
errors and the execution times obtained by our FMC procedure (top) and the 
Fluid Approximation of the mean behaviour (bottom) are reported in Tab. |T] 
Regarding the errors, we would like to remark that the Relative Errors (RE) of 
both the FMC and the Fluid Approximation reach their maximum in the very 
first instants of time, when the true satisfaction probability (i.e. the denominator 
of the REs) is indeed really small, but then they decay really fast as the values 
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Fig. 3. The satisfaction probability P(T) = Pf d (T) obtained by the Fluid Model 
Checking (left) and the Fluid Approximation of the mean behaviour (right) in the case 
N = 1000. The results are compared with those obtained by the DES. 


of P Fo ( t ) increase (this can be easily deduced from the values of the mean REs 
and the REs at final time). As expected, the accuracy of the approximations 
increases with the population size, and is already reasonably good for N = 500. 
Moreover, the resolution of the DDEs is computationally independent of N, and 
also much faster (approximatively 3 orders of magnitude in the case of the Fluid 
for N = 1000) than the simulation based method. Fig.[3]shows the results of the 
FMC and the Fluid Approximation in the case N=1000. 

5 Conclusions 


We defined a fast and efficient FMC procedure that accurately estimates the 
probability that a single agent inside a large collective system satisfies a time- 
bounded property specified by a single-clock DTA. The method requires the 
integration of a system of DDEs for the transient probability of an IMRP, and the 
exactness of the estimation is guaranteed in the limit of an infinite population. 


Future Work. During the experimental analysis, we realised that, on certain mod¬ 
els and properties, the DDEs (17 1 can be stiff, and their numerical integration in 
MATLAB is unstable (see also 8|). In the future, we want to address this issue 
by considering alternative integration methods 21 , investigating also numerical 
techniques for MRP with time-dependent rates 26 . Furthermore, we plan to 


prove the convergence of the Fluid Approximation of Sec. |3.1[ and to investigate 
higher-order estimates as in 


10 11 


Finally, we want to extend the FMC pro¬ 
cedure of this paper to validate requirements specified in the logic CSL Tj4 17 


and DTA properties endowed with multiple clocks (possibly considering the ap¬ 
proximation techniques defined in 19 and [4]). 
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A Proofs 


A.l Prob Z (jv> and Probz, and Measurability of X 




Consider a single agent of class si = ( S,E) in a population model X^ = 
(si and a timed property D = D (T) = (T, , I 5 , CC, Q, q 0 , F, ►). 

Let P 
in Sec. 


2.1 


t be the set of time-bounded paths of si accepted by D introduced 
and let Z' N >(t) and Z(t) be the two stochastic processes defined for 


the Fast Simulation in Sec. [2] 

The process Z^ N \t) is non-Markovian, and it becomes a CTMC only when 
considered in couple with the state X^ N \t) of X^ N \ Indeed, given the Markov 
process Y = (Y^ N \t), ..., Yjf\t)), that enlists all the stochastic processes that 
describe the state of the N agents in X^ N \ Z^ N \t) is defined to be the projection 
of Y on the component Y^ N \t) that represents the single agent that we are 
considering in our FMC procedure. For this reason, if we assume w.l.o.g. that 
i = 1, i.e. the stochastic process for the single agent is the first component of Y, 
the transient probabilities of Z^ N \t) and (Y^ N \t),... ,Y^\t)) are such that 


= k}= J2 P{(Y 1 W (t),...,Yi N \t)) = (k, S )}. 


( 10 ) 


ses "- 1 


Then, exploiting the equality (10), the probability measure Prob Z (N ) over path 
of ZW can be easily derived form Prob^ Z (N) X («)), which is the probability 
measure of the CTMC (Z^ N \ X^) defined in the standard way over cylinder 
sets of paths (cf e.g. [13]). 

The probability measure Prob z over the paths of the stochastic process Z(t ), 
instead, is the standard probability measure for ICTMC (cf. e.g. ID- 

Then, the measurability of d.t for Prob Z ( N ) and Prob z is just an adap¬ 
tation of Theorem 3.2 of [l3 to Z^\t) and Z(t). 


A.2 Convergence of the Approximation 

Consider a single agent of class si = (S,E ) in a population model X^ = 
(si,T^ N \x o^), an d a timed property D = O(T) = (T, , P$, CC, Q, qo, F, —»). 
Let Uj rf.n.T be the set of time-bounded paths of si accepted by D introduced 
and let Z^ N \t) and Z(t) be the two stochastic processes defined for 


in Sec. 


2.1 


the Fast Simulation in Sec. [2] Consider the probability measures Prob Z (N > and 
Prob z of Z( N \t) and Z(t), respectively, and define the satisfaction probabilities 
p( N )(T) = Prob Z (N) {Yjj,o ,t} and P(T) = Prob z Then the following 

theorem holds true. 

Theorem [3] For any T < +oo, 


lim P (JV) (T) = P(T). 

JV-H-oo 
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Proof. By a standard coupling argument, we can assume that both Z^ N ^ and Z 
are defined on the same sample space 12. Moreover, for each w £ 17, let Z^^lo, t ) 
and Z(oj,t) denote the state reached at time t in the path corresponding to u> 
of ZW(i) and Z(t), respectively. Then, according to Theorem [ 2 } for each time 

horizon T < + 00 , there exist a sequence cn £ M + , £n 0, such that 

Probes £ n | Vi < T, Z (N \u,t) = Z{u,t)} >l-e N . 

Now, fix T and N, let Et be the set of paths of total duration T for an agent 
of class si, and consider the (measurable) function xb : Et —> {0,1}, whose 
value xb(o') is 1 if the path a £ Et is accepted by D (i.e. the labels of er define a 
path of D that ends in one of the final states F), and 0 otherwise. By definition, 
p( N )(T) = E[xd(^^)] and P(T ) = E[xd(Z)], where Xb is evaluated on the 
paths of Z^ N } and Z restricted up to time T. Moreover, define 

= {w £ Q | zW{u,t) = Z(u,t), Vt £ [0,T]| , 


and l?o = 12 \ l?i. Then, xb (Z ( ' N ' ) ) = xb(Z) on 17 1 and Prob(f2 0 ) < ejv- Hence, 
we have 


E 


Xb(Z {N) )] -e[xb (Z) 


< E 


X0 (ZW) - xb(Z) 

Xd(Z^)- X d(Z) 




Xb 


(Z^)-xb(Z) 


dun ■ 


d^n < e-N - 0 . 


A.3 Results on the Product sin and the IMRP Z^ 0 (t) 

Consider a single agent of class si = {S,E) in a population model X^ = 
and a timed property D = D(T) = (T, , P$, CC, Q, go, F, —»). 

Let Ej j.n.T be the set of time-bounded paths of si accepted by D introduced 
and let Z^ N \t) and Z(t) be the two stochastic processes defined for 


2.1 


in Sec. 

the Fast Simulation in Sec. [2] Moreover, consider the product si n of agent class 
si and property D, and the IMRP Z^ B (t), both introduced in Sec. [3j 

Let AccPath(sio,T) be the set of time bounded paths of total duration T 
that are accepted by si b (i.e. such that they terminate in one of the final states 
of Mj)- Then, the result that guarantees that E^ t b,t is in 1:1 correspondence 
with AccPath(s/o,T) is just an adaptation to L^b, t and AccPath(s/o,T) of 
Lemma 3.9 in |l3|. 

Finally, the fact that P(T) = Probz W£ {AccPath[sii j,T)} comes from a (non¬ 
trivial) adaptation to time-dependent rates of Theorems 3.10 and 4.3 also in [13;. 
We plan to provide the details of this formal proof in the future journal version 
of this paper. 
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